{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "4ed546c8-cbad-4721-ba6b-7d185222990e",
   "metadata": {},
   "source": [
    "# Example 1\n",
    "No source, Single-connected domain, Dirichlet problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "228a690d-545c-48b5-a2d5-6a3e40d7963b",
   "metadata": {},
   "outputs": [],
   "source": [
    "import calfem.geometry as cfg\n",
    "import calfem.mesh as cfm\n",
    "import calfem.vis as cfv\n",
    "import numpy as np\n",
    "from toolkits import gauss_line_int\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b47a234-100f-4115-a5ee-670e2b590fdd",
   "metadata": {},
   "source": [
    "## 1. Geometry and Meshing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "1cacc716-7c70-4662-8bca-79ca9c312b4d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Info    : GMSH -> Python-module\n"
     ]
    }
   ],
   "source": [
    "g = cfg.Geometry()\n",
    "\n",
    "g.point([0.0, 0.0]) # point 0\n",
    "g.point([1.0, 0.0]) # point 1\n",
    "g.point([1.0, 1.0]) # point 2\n",
    "g.point([0.0, 1.0]) # point 3\n",
    "\n",
    "# g.spline([0, 1]) # line 0\n",
    "# g.spline([1, 2]) # line 1\n",
    "# g.spline([2, 3]) # line 2\n",
    "# g.spline([3, 0]) # line 3\n",
    "# g.surface([0, 1, 2, 3])\n",
    "\n",
    "factor = 20\n",
    "dir_id_1 = 100\n",
    "dir_id_2 = 200\n",
    "dir_id_3 = 300\n",
    "dir_id_4 = 400\n",
    "\n",
    "g.spline([0, 1], el_on_curve=factor, marker=dir_id_1) # line 0\n",
    "g.spline([1, 2], el_on_curve=factor, marker=dir_id_2) # line 1\n",
    "g.spline([2, 3], el_on_curve=factor, marker=dir_id_3) # line 2\n",
    "g.spline([3, 0], el_on_curve=factor, marker=dir_id_4) # line 3\n",
    "g.structuredSurface([0, 1, 2, 3])\n",
    "\n",
    "cfv.drawGeometry(g)\n",
    "cfv.showAndWait()\n",
    "\n",
    "mesh = cfm.GmshMesh(g,return_boundary_elements=True)\n",
    "\n",
    "mesh.elType = 2          # Degrees of freedom per node.\n",
    "mesh.dofsPerNode = 1     # Factor that changes element sizes.\n",
    "# mesh.elSizeFactor = 0.5 # Element size Factor\n",
    "\n",
    "coords, edof, dofs, bdofs, elementmarkers, boundaryElements  = mesh.create()\n",
    "\n",
    "element = edof - 1\n",
    "node = coords\n",
    "\n",
    "E = len(element)\n",
    "N = len(node)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "edf91826-efaf-416e-9306-d783aa7ea2c6",
   "metadata": {},
   "source": [
    "## 2. Boundary condition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "2ba4fdfb-bf49-49c4-b787-3c5e9f94e771",
   "metadata": {},
   "outputs": [],
   "source": [
    "dir_edge_1 = np.array([edge['node-number-list'] for edge in boundaryElements[dir_id_1]])-1\n",
    "dir_edge_2 = np.array([edge['node-number-list'] for edge in boundaryElements[dir_id_2]])-1\n",
    "dir_edge_3 = np.array([edge['node-number-list'] for edge in boundaryElements[dir_id_3]])-1\n",
    "dir_edge_4 = np.array([edge['node-number-list'] for edge in boundaryElements[dir_id_4]])-1\n",
    "\n",
    "# boundary edge with node in counterclock-wise order\n",
    "boundary = [edge['node-number-list'] for marker in boundaryElements for edge in boundaryElements[marker]]\n",
    "\n",
    "# def dir1(x, y):\n",
    "#     return x*(1-x)\n",
    "\n",
    "# def dir2(x, y):\n",
    "#     return -y*(1-y)\n",
    "\n",
    "# def dir3(x, y):\n",
    "#     return x*(1-x)\n",
    "\n",
    "# def dir4(x, y):\n",
    "#     return -y*(1-y)\n",
    "\n",
    "def dir1(x, y):\n",
    "    return np.zeros_like(x)\n",
    "\n",
    "def dir2(x, y):\n",
    "    return np.sinh(y)/np.sinh(1)\n",
    "\n",
    "def dir3(x, y):\n",
    "    return np.sin(x)/np.sin(1)\n",
    "\n",
    "def dir4(x, y):\n",
    "    return np.zeros_like(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d5fe1912-f8d3-4895-9e96-3156355eb71c",
   "metadata": {},
   "source": [
    "## 3. Assemble flexibility matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "fb18dbeb-89fc-4566-b9bc-7627f4bef33f",
   "metadata": {},
   "outputs": [],
   "source": [
    "F = np.zeros((N,N))\n",
    "\n",
    "for i, e in enumerate(element):\n",
    "    x1,y1 = node[e][0]\n",
    "    x2,y2 = node[e][1]\n",
    "    x3,y3 = node[e][2]\n",
    "    \n",
    "    A_e = 1/2*np.linalg.det(np.c_[node[e], np.ones((3,1))])\n",
    "    \n",
    "    L1 = np.sqrt((x2-x1)**2+(y2-y1)**2)\n",
    "    L2 = np.sqrt((x3-x2)**2+(y3-y2)**2)\n",
    "    L3 = np.sqrt((x1-x3)**2+(y1-y3)**2)\n",
    "    \n",
    "    left = np.linalg.inv([[(y2-y1)/L1,(x1-x2)/L1],\n",
    "                          [(y3-y2)/L2,(x2-x3)/L2]])\n",
    "    \n",
    "    Phi_L = np.c_[left, np.zeros((2,1))]\n",
    "    Phi_e = Phi_L  @ np.array([[-1/L1,1/L1,0],\n",
    "                               [0,-1/L2,1/L2],\n",
    "                               [1/L3,0,-1/L3]])\n",
    "    \n",
    "    Phi_e = (1/(2*A_e))*np.array([[x2-x3,x3-x1,x1-x2],\n",
    "                                  [y2-y3,y3-y1,y1-y2]])\n",
    "    \n",
    "    # print(np.linalg.norm(Phi_e-_Phi_e))\n",
    "    \n",
    "    F_e = A_e * Phi_e.T @ Phi_e\n",
    "    \n",
    "    F[np.ix_(e,e)] = F[np.ix_(e,e)] + F_e\n",
    "    # with np.printoptions(precision=2, suppress=True):\n",
    "    #     print(e)\n",
    "    #     print(F_e)\n",
    "    #     print(F)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7a01a3f7-c7e0-40b2-aeef-6ac3d2787a52",
   "metadata": {},
   "source": [
    "## 4. Assemble kinematic vector"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "c855dbeb-fe00-4b49-8f23-8514e2fd3e53",
   "metadata": {},
   "outputs": [],
   "source": [
    "du = np.zeros(N)\n",
    "\n",
    "for edge in dir_edge_1:\n",
    "    i, j = edge\n",
    "    L = np.linalg.norm(node[i]-node[j])\n",
    "    up = gauss_line_int(dir1, *node[i], *node[j])\n",
    "    du[i] = du[i]-up/L\n",
    "    du[j] = du[j]+up/L\n",
    "\n",
    "for edge in dir_edge_2:\n",
    "    i, j = edge\n",
    "    L = np.linalg.norm(node[i]-node[j])\n",
    "    up = gauss_line_int(dir2, *node[i], *node[j])\n",
    "    du[i] = du[i]-up/L\n",
    "    du[j] = du[j]+up/L\n",
    "\n",
    "for edge in dir_edge_3:\n",
    "    i, j = edge\n",
    "    L = np.linalg.norm(node[i]-node[j])\n",
    "    up = gauss_line_int(dir3, *node[i], *node[j])\n",
    "    du[i] = du[i]-up/L\n",
    "    du[j] = du[j]+up/L\n",
    "\n",
    "for edge in dir_edge_4:\n",
    "    i, j = edge\n",
    "    L = np.linalg.norm(node[i]-node[j])\n",
    "    up = gauss_line_int(dir4, *node[i], *node[j])\n",
    "    du[i] = du[i]-up/L\n",
    "    du[j] = du[j]+up/L\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5619f835-d960-4905-a43c-5b77ff0dcc40",
   "metadata": {},
   "source": [
    "## 5. Constraint"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "0dba5f60-fa70-4d41-ac50-72d386d3aa69",
   "metadata": {},
   "outputs": [],
   "source": [
    "FF = F[:]\n",
    "ddu = du[:]\n",
    "FF[-1,:] = 0\n",
    "FF[:,-1] = 0\n",
    "FF[-1,-1] = 1\n",
    "ddu[-1] = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4c4ffc26-a9f0-408c-80eb-723329c3b8ee",
   "metadata": {},
   "source": [
    "## 6. Solve and show"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "975aa11a-a4b7-49fb-869b-39843c1242e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "s = np.linalg.inv(FF) @ ddu"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8e4ab57-f7ab-4b05-a27f-61405d7d3079",
   "metadata": {},
   "source": [
    "## 7. Construct the flux field"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "d4a4a3c5-081a-43ff-9f2a-1603428996bf",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<visvis.wibjects.colorWibjects.Colorbar object at 0x0000021FAF8B6280>\n",
      "ERROR calling '_OnAxesPositionChange':\n",
      "  File \"D:\\20220326(EFEM-Poisson)\\EFEM\\EFEM\\lib\\site-packages\\visvis\\wibjects\\colorWibjects.py\", line 622, in _OnAxesPositionChange\n",
      "    self.position.x = axes.position.width+5\n",
      "  File \"D:\\20220326(EFEM-Poisson)\\EFEM\\EFEM\\lib\\site-packages\\visvis\\core\\misc.py\", line 217, in fsetWithDraw\n",
      "    fset(self, *args)\n",
      "  File \"D:\\20220326(EFEM-Poisson)\\EFEM\\EFEM\\lib\\site-packages\\visvis\\core\\base.py\", line 1126, in fset\n",
      "    self._Update()\n",
      "  File \"D:\\20220326(EFEM-Poisson)\\EFEM\\EFEM\\lib\\site-packages\\visvis\\core\\base.py\", line 928, in _Update\n",
      "    self._CalculateInPixels()\n",
      "  File \"D:\\20220326(EFEM-Poisson)\\EFEM\\EFEM\\lib\\site-packages\\visvis\\core\\base.py\", line 997, in _CalculateInPixels\n",
      "    raise Exception(\"Can only calculate the position in pixels\"+\n",
      "Exception: Can only calculate the position in pixels if the owner has a parent!\n",
      "\n"
     ]
    }
   ],
   "source": [
    "q_h = np.zeros((E,2)) \n",
    "\n",
    "for i, e in enumerate(element):\n",
    "    x1,y1 = node[e][0]\n",
    "    x2,y2 = node[e][1]\n",
    "    x3,y3 = node[e][2]\n",
    "    \n",
    "    A_e = 1/2*np.linalg.det(np.c_[node[e], np.ones((3,1))])\n",
    "    \n",
    "    Phi_e = 1/(2*A_e)*np.array([[x2-x3,x3-x1,x1-x2],\n",
    "                                [y2-y3,y3-y1,y1-y2]])\n",
    "    \n",
    "    q_e = Phi_e @ s[e]\n",
    "    q_h[i] = q_e\n",
    "\n",
    "cfv.figure()\n",
    "\n",
    "# Draw the mesh.\n",
    "\n",
    "cfv.drawElementValues(\n",
    "    ev=q_h[:,1],\n",
    "    coords=coords,\n",
    "    edof=edof,\n",
    "    dofs_per_node=mesh.dofsPerNode,\n",
    "    el_type=mesh.elType,\n",
    "    title=\"Example 01\"\n",
    "        )\n",
    "cfv.showAndWait()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d79ef90-62c8-4a33-9848-c3d5ce5b4d9b",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "EFEM",
   "language": "python",
   "name": "efem"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
